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O ■ Abstract 

o 

C3 ' Lectures on Lattice Field Theory and Lattice QCD given at the 

Graduate Students Association (GSA) Summer School (Fermilab). 
In these lectures we provide a short introduction to the Monte 
^ i' Carlo integration method and its applications. We show how the origin 

^ • of ultraviolet divergences if Field Theories is in the undefined formal 

product of distributions and how one can define the Path Integral in 
terms of regularized distributions in order to cancel these divergences. 
^ ' This technique provides the only non perturbative regularization pro- 

cedure of continuum Field Theories and, at the same time, provides 
a practical method to compute correlation (Green) functions (using 
Monte Carlo integration for the regularized path integrals). We then 
apply these tools to formulate QCD on a lattice. Some of the exam- 
ples are accompanied by complete computer programs. 

Freely download libraries and examples from: 

http : //latticeqcd . f nal . gov/sof tware/f eriiiiqcd/| 



Introduction 

Two are the main tasks a physicist has to confront with: 

• [Induction] : given the symmetries of the measured observables, build 
the underlying theory (i.e. write down an action). 

• [Deduction]: given the action, S, compute correlation functions and, 
from them, physical observables (to test the theory and to make pre- 
dictions). 

In this notes we will focus the second task. 

In the first section we wil see how the perturbative expansion of the Path 
Integral does not provide a satisfactory definition of the latter and a non- 
pertubartive regularization is necessary in order to define it properly. 

For us the word "non-perturbative" means "exact up to a given precison 
that can be arbitrarily small". 

In the second section we will show how to compute numerically K- 
dimensional integrals (for large integer K) using Monte Carlo techniques. 

In the third section we will define the regularized Path Integral in terms 
of regularized distributions. This is equivalent to discretize the the space- 
time on which the quantum fields are defined. The lattice spacing a will play 
the role of an ultraviolet cut-off. We will then define the continuum Path 
Integral in terms of i^- dimensional integrals in the limit a ^ 0. We will 
show how this limit is the origin of ultraviolet divergences and how one can 
renormalize the theory by giving an a dependence to the coupling constants 
that appear in the action. 

Finally in the third section we will see how one can approximate a con- 
tinuum Path Integral of QCD with a i^-dimensional integral (for a finite K) 
and compute it numerically using Monte Carlo integration. We will discuss 
the sources of numerical errors and we will present, as an example, one full 
Lattice QCD application (the computation of fB^/iT^B)- 

Some more examples of typical Lattice QCD computations are given in 
the Appendix. 

The emphasis in these lectures will be given to three aspects that we 
consider crucial and make of lattice a privileged tool in respect to other 
model independent methods: 

• The lattice regularization provides the only non-perturbative definition 
of Path Integral and, therefore, of Field Theories. 



m 



• It is possible to quantify with precision the error committed in the 
numerical approximation of the integrals. 

• It is, in principle, possible to reduce arbitrarily this error by approach- 
ing the continuum limit (reducing the lattice spacing) and increasing 
the statistical samples in the Monte Carlo integration. 

This lectures are intended to be an introductory tutorial and they are 
not meant to be complete and/or exhaustive on the subject of lattice QCD. 
For more complete introductory reviews on the subject see |]1| and 0. For a 
more complete and formal review see 0. 

I wish to thank G. Chiodini, B. Dobrescu, E. Eichten, J. Juge, A. Kro- 
nfeld, P. Mackenzie and J. Simone for helpful comments and suggestions 
regarding these notes. 

This work was performed at Fermilab, a U.S. Department of Energy Lab 
(operated by the University Research Association, Inc.), under contract DE- 
AC02-76CHO3000. 
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1 Correlation functions, masses and matrix 
elements 



He who loves practice without theory is 
like the sailor who boards ship without 
a rudder and compass and never knows 
where he may cast 

(Leonardo da Vinci) 



The symbol of Path Integral |^ 

(O|r{0(xi)...0(x„)}|O) =^ / [d0]0(xi)...0(xOe-'^-['^J (1) 



provides a definition of the most general Euclidean correlation function (the 
left hand side) in terms of an infinite-dimensional integral (the path integral 
at the right hand side). 5e[0] is the Euclidean action of the system and (j){x) 
represents the degrees of freedom of the system as function of the space-time 
Euclidean coordinates. 

Without loss of generality we will only deal with Path Integrals in the 
Euclidean space since Minkowskian correlation functions can be obtained by 
analytic continuation of the Euclidean ones. In particular we are interested 
in extracting particle masses and matrix elements which do not have an 
explicit time dependence. Therefore, up to a phase, they are the same in the 
Minkowskian and in the Euclidean formulation of the theory. 

In the case of the 0^ scalar model the perturbative expansion, after renor- 
malization, is always convergent therefore it provides us with a definition of 
the Path Integral. Unfortunately... 

In general, the perturbation expansion does not provide a satis- 
factory definition of the Path Integral. 

Let's consider, for example, the case of QCD. The perturbative expansion 
of Path Integrals of QCD is only convergent at high energy {p >> 200MeV), 
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i.e. at short distance (a = 1/p « Ifm ). In fact the expansion parameter 
Us = ^^^ is scale dependent and it becomes big at low energy, i.e. large 
distance |^, §], • Moreover usual regularization prescriptions (such as Pauli- 
Villars and dimensional regularization), which are necessary to give sense to 
the Path Integral, are only defined in a perturbative way. 

One could argue that, even if the perturbative expansion of eq.(|l|) is 
only valid at high energy, it still defines eq.(P because one can analytically 
continue the correlation functions from high to low energy and make somehow 
sense out of them. This is not correct! 

First, the perturbative series is an "asymptotic series" and it does not 
necessarily converge. Second (the most important point), the perturbative 
expansion is performed around the classical minimum of the action. The Eu- 
clidean action of QCD has many different minima that correspond to multi- 
instanton configurations. These are tunneling transitions between different 
vacua of QCD. These effects are not taken into account in the perturbative 
expansion but, nevertheless, they play a crucial role at low energy. They are 
believed to be responsible for the phenomenon of confinement IQ . 

The lattice discretization of the space-time provides us with the 
only non-perturbative regularization of the Path Integral and, there- 
fore, with an EXACT definition of the latter. 

For a small and finite cut-off (the lattice spacing a) the regular- 
ized theory can he seen as an effective theory of the continuum 
one (corresponding to the limit a ^ Q). Such an effective the- 
ory is finite and free of ultraviolet divergences therefore it can he 
simulated numerically! 

Moreover, in the case of QCD and gauge theories in general, lat- 
tice regularization has the nice feature to preserve gauge invari- 
ance even for every finite a. 

We now introduce two examples in which we require a non-perturbative 
definition of the Path Integral to compute some important phenomenological 
quantities. 



1.1 2-point correlation functions 

QCD is the theory of strong interactions, hence it must predict the masses of 
mesons and baryons from first principles. Let's consider here, for example, a 
B meson. 

The following current 

J^'ix) = h{x)Yl^q{x) (2) 

(where x^ = (xq = tx,xi,X2,X3), tx is the Wick rotated time and h (q) is 
the fermionic field representing the heavy (light) quark) has the same flavor 
quantum numbers of a B meson. If we apply J° to the vacuum, it must 
create some linear combination of a static B meson and its excited states 

J°|0) = eo\B) + £i|5«) + S2\B^^^) + ... (3) 

(in general there is a continuum of states but, for notation, we write them here 
as a sum.) From now on all our states are normalized as {B^"-^\B^^^) = 2mg(n) . 
We define the following (zero momentum Fourier transform of the) two 
point correlation function 



C^itx) = / d3x(0|J°(x)j0t(0)|0) 

[dA^] [dq,] [dgj fj d'^j\x)j'\0)] e-5-^'° (4) 

One can insert in the correlator a complete set of states... 

d-^x(0|J°(a:)^-^^ ^J°^(0)|0) 

2m5(„) 

\ ]3{n)\ / f){n)\ 

= ^(0|J°(0)e-^*- '% ^^^ ' j°t(o)|o) 

n 

^ ^ |(0|J°(0)|i?("))P ^_^(„H^ 



2mi 



n 



= \Zo\'e-^"'^ + \Zi\'e-^"'^ + ... (5) 

where e""^*"" is the Euclidean translation operator and 

if|5(")) = ^(")|5(")) (6) 



The long distance behavior of C2{tx) is dominated by the exponential 
associated with the lightest state |i?*^°^) (with energy E^^^ = niB) 

C^iQ ~ iZol^e-"^^*^ (7) 



Therefore by definition^] 

^ _ (0|J°(0)|i?W) _ 1 



fBy/rn^ (8) 



Eq.(^) is not defined (yet) in the limit tx ^ oo because this limit is clearly 
large distance and therefore non-perturbative. 

If we had a non perturhative definition of eq. (^, we could compute 
C2{tx) o,nd extract both itlb and fs by fitting the result with eq. (0j. 

An example of two-point correlation function computed on a lattice (using 
program C2 . C in the Appendix) is shown in fig. |]. 

1.2 3-point correlation functions 

Vcb is a very important parameter of the Standard ModelQ. It must be de- 
termined by comparing QCD contributions with experiment. By definition 
i0: 



\Vcb\^ = -, — -= ^-^ — ^ [perturbative factor] [experiment] (9) 

From C2{tx) we have the coefficients in the numerator but we also need a 
QCD prediction for the matrix element in the denominator. To reach this 
goal we define a three-point correlation function 

C^{tx.ty) = Jd^^Jd'y{0\J%x)OJ%-y)\0) 

[dA,][dq,][dq,] (^j d'^ j d'yJ\x)OJ\-y)\ e-^HlO) 



1 {Q\J^^\B)=fBP^' 

"^Vcb is one of the elements of the Cabibbo-Kobayashi-Maskawa matrix V and it is 
known with very poor precision. A non unitary CKM matrix would be a clear signal of 
physics beyond the Standard Model. 
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Figure 1: Example of a numerical two-point correlation function for a B 
meson for some (arbitrary) input values of the light and heavy quark masses. 
The plot is symetric arout t^. = 8 because of periodic boundary conditions. 
The curve (for t^ < 8) can be fitted with C2{t) = \f'^mBe~™'^^ to determine 
TUB and JB- 



where 

O = h{0)YLq{0)h{0)YLq{0) (11) 

(L, defined in the appendix, is the left-handed projector) We play the same 
trick as before (inserting two complete sets of states) and we obtain 

^-^ J J Zm^in) Zm^im) 

n,m " " Lj Lj 

n,m V u u 

which for t^, -^ oo is dominated by 

Cit^^ty) ^ |^^|2,-r..(t.+*.)M^M (13) 

'^ ta:,ty-*OC 2niB 

Again... 

// we had a non-perturbative definition of eq. ^T^), we could com- 
pute C^itx.ty) and extract matrix elements (such as {B\0\B)) by 
fitting the result with eq.^Jl 



We will see in the next sections how the lattice provides both a definition 
and a practical way of computing the right-hand side of eq.(|l]), (H) and (plQl). 



■'This is what one actually does to determine Vet- At present the main error on Vcb is 
due to theoretical uncertanties. 



2 Monte Carlo Integration 



Although this may seem a paradox, all 
exact science is dominated by the idea 
of approximation 

(Bertrand Russell) 



In this section it will be shown how to compute numerically multidimen- 
sional integrals using a statistical method (the Monte Carlo integration). 

2.1 Riemann integrable functions 

Let's consider the simplest of the integrals 



f{x)dx (14) 

Riemann gave us a definition of this integral: 

To check if the function f{x) in the domain V = [a,f3] is Riemann inte- 
grable we divide the domain in N small intervals 

[xq = a, xi], [a;i, 0:2], ..., [xn-2, xn~i], [xn~i,xn = /?] (15) 

and compute 

N-l N-l 

£{N) = a 2, max f{x) — a >^ min f{x) (16) 



. X^\Xi,Xi^]\ ■^— ^ X^\Xi,XiJ^x\ 

1=0 1=0 



where a =^ ^ 



/(t) is Riemann integrable if e{N) goes to zero for N -^ 00, i.e. 
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Figure 2: Convergent behavior of different numerical discretizations for the 
integral of eq. (|T^ . 



This definition of integration is also very practical because it gives a 
method of performing numerical integrations: 






There are other possible ways of approximating a continuum integral with 
a discrete sumQ and, if they converge, they all converge to the same number. 

The difference between the different integration methods resides 
in the behavior of the numerical integral as function of N: some 
methods converge faster than others. For practical purposes this 
is a crucial issue. 

programl . c is a C program that computes the integral 

I = sin(7ra;)da; ~ 0. 63662 (19) 

Jo 

using three different simple methods. The convergent behavior of the differ- 
ent methods is shown in fig.l. The error done when one stops at a finite A^ 
is known as disctretization error. 

For multidimensional integrals this method turns out to be too slow for 
practical applications. Therefore we need a new tool, Monte Carlo integra- 
tion. 

2.2 Basic Monte Carlo 

Monte Carlo integration is another numerical method for integration and it 
has a statistical foundation. The algorithm is the following: 

^for example using trapezoids instead of rectangles (Newton-Cotes method) 

[^ f{x)dx = X;' « • ^(^')+/^'''^'^ + Oif'/N') (17) 

•^° i=0 ^ 

or using interpolating polynomia (Simpson's method) 

/' /(x)d. = X: a • /(^-)+4/(^'+i) + /(^^+2) ^ o{r/N^) (18) 

See the Numerical Recipes for more details. 



• Generate a set of A^ random points {x'*'} with uniform distribution in 
the integration domain T> 

• For each point xf*^ compute /(x'^J) 

• Compute the average I{N) = ^^ X]j=o /(^^'^O 

If f{x) is Riemann integrable, the function I{N) converges to the 
integral of f{x) when N ^ oo. 

program2 . c is a C program that computes numerically, using Monte Carlo 
integration, the same integral of eq. (p^Q]). Its convergent behavior is shown 
in fig.2 and compared with the standard methods showed in section 2.1. 

This algorithm can easily be extended to arbitrary dimensions. For ex- 
ample programs . c computes the following integral 

/ = / dxo / dxi / dx2{Sxlxi + 2x1) — ^ (20) 

Jo Jo Jo 

and its convergent behavior is shown in fig.3. Monte Carlo integration is not 
very efficient for ID integrals but it becomes more and more efficient (when 
compared with conventional integration methods) for higher dimensional in- 
tegrals. 

From now on we will use the upper index x'*' to label the Monte Carlo 
points and a lower index Xj to label a particular coordinate of a given point. 

[31 

For example Xg is the second coordinate of the 3rd point. 

2.3 Metropolis Monte Carlo 

In the next section we will deal with f^-dimensional integrals of the form 

4 = / dxo / dxi... / (ixi4'_i/(x)P(x) (21) 

Jo Jo Jo 

where x = (xq, xi, ..., xk-i) (22) 

and we have to compute them for different functions /(x) but the same 
P(x) 0. The clever trick we play is to modify our algorithm in the following 

way: 

^think for example to the case of many Euclidean green functions Gk defined as 

Gfc= /"[dx]/fc(x)P(x) (23) 

10 



Monte Carlo Integration Method 
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Figure 3: Convergent behavior of Monte Carlo Integration of eq.dl 
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Multidimensional Monte Carlo Integration 
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Figure 4: Convergent behavior Monte Carlo Integration of eg . (pO|) . 
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Generate a set of A^ random points {xW} in the integration domain 
V using P(xW) for the probabihty distribution of the point x^ = 

For each point x^ compute /(x^) 
Compute the average 

7V-1 



j=0 

where Vol(T') is the volume of the integration domain. 

// /(x)P(x) is Riemann integrable then the function I{N) con- 
verges to the integral eq.^\) when N —* oo. The factor -P(x) in 
the integrand has been absorbed into the probability distribution 
of generating the random points. 

We are now left with the problem of generating random points x^ with a 
given distribution P(xW). There are many algorithms that do this job. The 



simplest one is the Metropolis algorithm [|TT 



1. Start with i = and a point x^ chosen at random in the integration 
domain. 

2. Generate another random point y in the integration domain and a 
random number a in the interval [0, 1). 

3. If P(y)/P(xH) > a then x'^+^l = y else x'^+^l = x'^l. 

4. Increase i hj 1 and repeat steps 2,3,4. 

The succession of points {x^} generated by this algorithm has the re- 
quired probability P(xW). This kind of succession is also known as Markov 
chain. 

Each point of the Markov chain is called a configuration. 



where P(x) ex e ■^s^'') and Se{x) is the Euchdean action of a system in the configuration 

X. 
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Figure 5: Convergent behavior of the Metropohs algorithm. 
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The Metropolis Monte Carlo only works with a real probability function, 
P{x), but the functions f{x), and the integration variables x, can be complex, 
vectors or tensors. 

program4 . c computes the following 3D integral^ 

/ = / dxo dxi I rfx2/(x)P(x) (25) 

Jq Jo Jo 



where P(x) = e-("o+-?+-i) 
and /(x) = 12^x1x1x2 

Its convergent behavior is shown in fig.4. 

2.4 Bootstrap errors 

Monte Carlo integration does not just provide us with a way to compute 
multidimensional integrals. There are two techniques, known and Jackknife 
and Bootstrap [|T2|, that permit us to evaluate the statistical error we commit 
when we truncate the Markov chain and we compute the numerical integral 
on a finite number of configurations, A^. We consider here the Bootstrap 
method. 

Let's assume we have a finite number, N , of configurations that constitute 
the Markov chain {x^} and the integral / which is computed numerically as 

/^/(iV)t/M^^V(x[^') (26) 

i=0 

Since in general the /(x^) are not Gaussian distributed, standard 
error analysis is not applicable. 

The Bootstrap algorithm consists of the following steps: 

• Construct a table of A^ x M integer random numbers kij where i G 
{0,l,...,Ar-l}, j e {0,l,...,M-l}andA:ij G {0, 1, ..., A^- 1} for each 
couple (i, j) (M is an input parameter that we choose to be equal to 
100). 



^it can easily be computed analytically by factorization and one finds / = 2.43131 
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For each j compute 

N-l 



i=0 

• Reorder the set {Ij} so that, for each j, Ij < Jj+i 

The result for the integral I lies between 1^3 and Iqq at 65% con- 
fidence level. 

The idea behind the Bootstrap algorithm is that, for each j, and for 
N » 1, the j-th Markov chain {x''^''^]} (chain in i) has the same prob- 
abihty distribution of the original chain {x^}. Therefore the different Ij 
(the numerical integral computed using the derived chain j) have the same 
probability distribution as that of the derived chain j associated to Ij. 

programs . c implements the Bootstrap algorithm as complement of the 
Metropolis algorithm. Fig. 5 shows the same data of fig.4 including the Boot- 
strap error. The kind of error one does when truncates the Markov chain is 
known as statistical error. 
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Figure 6: Convergent behavior of the Metropohs and Bootstrap algorithm. 
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3 Definition of Path Integrals 



Calculus required continuity, and con- 
tinuity was supposed to require the in- 
finitely little; but nobody could discover 
what the infinitely little might be 

(Bertrand Russell) 



The concepts of Regularization and Renormalization play a fundamental 
role in the definition of the Path Integral (and in physics in general) 
There is a physical reason for it: 

• One never measures the value of a field (associated to a particle) in 
every point in space-time but one measures its integral over the test 
function of the physical detectors, which have a finite extension. There- 
fore there are mathematical reasons to require that the fields are defined 
in the space of distributions. 

• One wants to model the unknown short distance physics by introducing 
a Lagrangian density which contains only local (contact) interactions. 

If one tries to combine the previous statements in a Quantum Field The- 
ory, one encounters the problem of divergences and must find a way of dealing 
with them. The mathematical origin of these divergences is the presence of 
undefined formal products of distributions in the Lagrangian from which the 
path integral is computed. Regularization and Renormalization are in fact, 
in mathematical terms, the solution to the problem of defining these products 
of distributions!]. 

We analyze here, as an explanatory example, the problem of defining 
the product of 6 functions by regularizing them by a sequence of smooth 
functions that become more localized at zero. We then show how an arbitrary 
quantum field can be expanded and regularized using delta functions and 



^For general reviews on the subject see |1^, Q. For an application to Quantum Me- 
chanics see Q 
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how the presence of a finite spatial cut-off in the regularized distributions is 
equivalent to a finite cut-off in the momentum expansion of the field itself. 

This is not intended to be an introduction to Renormalization, but it is 
presented as an alternative view of its meaning having in mind the lattice as 
typical regulator. 

3.1 Toy example: Regularizing distributions 

The 6 function is a distribution which is defined as 

/ 5{x - x)F{x)dx = F{x) (28) 

for any smooth test function F{x). The same S function can be thought of 
as the limit of an ordinary function 6{a, x) 

5{x) =\im5{a,x) (29) 

where 5{a, x) must be enough regular, localized in x within a precision a and 
its integral normalized to one. This procedure is called regularization. Some 
possible regularization schemes are0 

S{a,x) = -[e{x + a/2)-e{x-a/2)] (30) 

5{a,x) = (vra^) ^ exp(— x^/a^) (31) 

6ia,x) = ^^"("^/"^ (32) 

TTX 

They are sketched in figure ^ The different schemes are equivalent in the 
sense that they give the same result for the following limit 

lim / 6{a, x — x)F{x)dx = F{x) (33) 

a^O J 

but they do not give a well defined limit in expressions of the form 

lim [[6{a,x-x)r+^F{x)dx=7 (34) 

a^O J 

^The 6{x) function is defined to be for a; < and 1 for a; > 0. It is discontinuous in 
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Figure 7: Examples of possible regularizations for a delta function. The x 
axis is in units of a, the y axis is in units of a~^. 



In fact eq. (^) is divergent as a~". Suppose one wants to give a well defined 
meaning to the limit in eq. (|3^) by removing somehow the divergence that 
occurs. One way of doing it is by making F{x) dependent on a, according 
with the prescription for 5(a,x), 



F{x)^FR{a,x) 



Z-^-{a)F{x) 



(35) 



where Z^{a) = a -[- 0{a^). In other words the divergence of the integral is 
absorbed in the normalization of the function F{x). This procedure is called 
renormalization and it depends on which regularization has been chosen. 
From now on the scheme of eq. ( |30D will be considered in particular. After 
renor malizat ion 



/"[(5(a, X - x)]"+^Z^"(a)F(x)dx = const. + 0{c 



(36) 



and its limit for a ^ becomes well defined. Therefore, up to order a terms 
one can redefine the integral of eq. (53) in the following way 



[5{x - x)f^^F{x)dx =^lim 



[(5(a, X — x)] 



T-^NTn+l rv — n 



Z^'^{a)F{x)dx (37) 



The situation can be even more complicated if F{x) itself is defined in terms 
of delta functions. For example one can consider the case when F{x) = 
exp[gS'^{x)]. In this case it is not sufficient to regularize 6 and renormalize F 
to get rid of the divergence, one is forced to renormalize g as well. 

It is a general statement that, if the function F{x,g) depends on some 
constant g, one has to renormalize the constant 



g -^ 9R{a) 



(38) 
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by imposing a constraint 

/ [6{a, X - x)]"+^Z^"(a)F(x, gR{a))dx = const. (39) 

Eq. (|39|) fixes the behavior of gnia) as function of a. Its solution, gnia), 
can have a non-trivial behavior in a. Eq.(p9D is a particular case of what 
is generally known as the Renormalization Group Equation (RGE) [Q. The 
behavior oi gfi{a) versus a is called running. Another common way of writing 
the renormalization group equation is 



d 



dloga 
or explicitly 



[6{a, X - x)]"+^Z^''{a)F{x, gR{a))dx = (40) 



(a-^ - P{9r)^ + niign)] j i^a, x - x)r+'Z],^F{x, gR)dx = (41) 



where 



(3{9r) = - ^1 gnia) 

a log a 

a log a 



(42) 

9r 

(43) 



9r 

If the original constant g is dimensionless, gR{a) must also depend on some 
other scale, say A, to cancel the dimension of a. In other words gR must be 
a function of aA, an adimensional quantity. This simple example shows how 
the renormalization procedure may force one to introduce a second scale A 
of which the renormalized constant is a function. This phenomenon is called 
dimensional transmutation. 

Dimensional transmutation is also a characteristic of QCD (and gauge 
theories in general). In typical physical problems, a is a free parameter and 
it can be chosen (the physics does not depend on it providing it is small 
enough). A, on the other side, characterizes the typical scale of the physics 
one is describing. If one measures gR{aA) at some arbitrary physical scale 
a = a one can uniquely determine the value of A. Once A is known one can 
predict the value of gR at any other scale a. 
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Figure 8: [left]: Example of a possible regularizations for the delta function. 
The X axis is in units of a, the y axis is in units of a~^. [center-right]: Example 
of a continuum field configuration (p and its approximation with a linear 
combinations of regularized delta functions. 

3.2 Defining the Path Integral 

We now go back to the most general correlation function, defined in terms 
of the Path Integral. We rewrite it as 



(O|T{0(xi)...0(xO}|O) 



def 



[d4>]F[(P{x),g] 



where F[(f){x),g] is the integrand 

F[0(x),(7]t/0(a;O...0(x„)e-^-l<^'^] 



(44) 



(45) 



and g is the coupling constant that appears in the action. For the moment 
we simply consider a one-dimensional scalar field theory 0(x) defined in the 
interval [0,L]. The specific form of the action, SE[4>,g] is unimportant but 
we assume that the action contains an interaction term of the form 



g / 0"(x)dx 



(46) 



with n > 2. This makes F[(j),g] a non-trivial functional of 



X 



"Lattice regularization" is the way to regularize the integral ^4li ) 
by approximating the fields with sums of regularized distributions. 
This is equivalent to discretize the space-time on which the fields 
are defined. 
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For a = 0.1 (K = 10): 



[d0]F[0, g] t/ j d0o... y d09F[0, g] 



For a = 0.05 {K = 20): 

[d<P]F[<P,g] 



+ F[ 




+ ... + F[ 



F[<P,g] 




+ ... + F[ 



For a = 0.025 {K = 40): 

[d<P]F[<P,g] 



F[<P,g] 




+ ... + F[ 



And at the continuum limit, a -^ {K -^ oo) 

[d0]F[0, (?] =^ lim y d0o... I d^K-i F[0, ^] 




+ ... + F[ ;;; 



(47) 



(48) 



(49) 



(50) 



J (?i"(a::)da; is finite 



J (f)"{x)dx is finite 



V 

<p"{x)dx^oo 



Figure 9: Example of lattice regularization of the Path Integral 
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In practice one introduces a niimimum lenght scale, the "lattice spacing" 
a, and expresses the most general path (field configuration) 0(x) in terms of 
regularized delta functions (as shown in fig. pi). 



K-l 
def 



x) ^ (f)^^^\a, x) = ^ (j)kS{a, X - ka) (51) 



fc=0 

where 

-ka+a 



def 
>fc = 



- / (l){x)dx (52) 

0- Jka 



are ordinary integration variables and 5{a, x) is a regularized delta func- 
tion (it is normalized to one and peaked at zero within the length scale a, 
fig-ileft]). 

With this prescription, for every finite a = L/K, one can define 

"[d0]F[0(a;),^]t/yd0oy"d0i...y"d0^_iF[0i^"(a,x),^] + O(a) (53) 

^^ V ' 

K^L/a 

Note the upper index (a) which identifies the regularized Path Integral 
with lattice spacing set to a. 

Due to discretization, the continuum variable x has been replaced by the 
discrete index k. 

For every finite K, the right-hand side of eq. (^) can be computed using 
the Monte Carlo technique discussed in the preceding chapter. In this case 
the integration variables 0^ replace the Xk of the last chapter. 



The regularized Path Integral, eq. ^5^), is well defined for any 
finite a. The limit K -^ oo is infested by divergences and we 
must give a prescription to make sense to this limit. 

Divergences associated with limit a ^- Q^K ^> oo (at L = Ka = constant) 
are called ultraviolet, while those associated with limit K,L -^ oo 
(at a = L/K = constant) are called infrared. 

Fig. H shows in a schematic way, how the path integral, eq. (|50| ) can can 
be approximated by finite multidimensional integrals, eqs. (|47| - ^9D , and the 

24 



fields are defined on a lattice. For each finite a, the "infinite sum on the path" 
(eqs. p7| - ^9D ) is well defined since all the divergences that may appear can be 
absorbed in the normalization of the integration measure and, for each path 
0, the functional F[0, ^f] is finite. In the limit a — > the regularized delta 
function S{x, a) become more and more peaked and approaches a Dirac 6{x) 
function (eq. (|50|)). Since the integrand F[(f){x),g] is non-linear in the field 
and the product of delta functions is not defined, F[(f),g] diverges on those 
configurations 4>{x) ~ 6{x). Therefore the Path Integral diverges. 

These divergences are not physical since, in any practical experiment, 
one only has a finite resolution, a, and one cannot discriminate a S{x) form 
6{x, a) if a < a. This is why one is allowed to substract these divergences or 
ignore them by considering the regularized theory an "effective theory" . 

The way out is the following: 

To have a well defined limit a ^ 0,K —>■ oo (at L = Ka = constant) 
one must impose that the result of the regularized path integral is 
independent from a. The only way to do it is to make the field 
normalization and the coupling constant (the g of eq. ^jB^)) de- 
pendent on a 

g^gji{a,A) (54) 

(the constant A must be introduced because in general a and g do 
not have the same dimensions). This makes the physics indepen- 
dent by the lattice scale a. 

One does it by choosing a particular correlation function (identified by 
the functional integrand F[(j),g] and imposing the contraint 



d log a 



d0o / d0i... / d^K-iF[4>{x),gn{a,A)] 



K~L/a 



(55) 



This determines the behavior (the running) of gR{a). Eq. (^5|) is nothing 
else than a way to write the Renormalization Group Equation for a lattice 
regularized theory. The appearence of A is called dimensional trasmutation. 
The effects of very short distance physics (below the length scale a) only 
appear in the regularized correlation functions through corrections which are 
proportional to a, or through the renormalization of the coupling (and, in 
general, of the masses) |jT9|. 
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It should be remarked that this procedure of defining the limit 
a ^ cannot be carried out for an arbitrary theory since there 
may be more sources of divergences than coupling constants. If 
this limit can be defined the theory is said to be renormalizable. 

Usually one distinguishes between the "bare" parameters that appear 
in the regularized Lagrangian (for a finite value of the cut-off, a) and the 
"dressed", "physical" or "renormalized" parameters that are defined, i.e. 
measured, by actual experiments. If one takes the limit a — >^ 0, the bare pa- 
rameters lose any physical meaning and one must carefully define the renor- 
malized ones (one says to choose a prescription). If one is happy of keeping 
the cut-off small but finite one is allowed to identify the renormalized and the 
bare parameters, because these can now be measured. This is the approach 
one uses on the lattice and it corresponds to the Kadanoff- Wilson approach 
to renormalization. 

One can write a RGE both for the bare parameters (as function of the 
cut-off) or, equivalently, for the renormalized ones (as function of the renor- 
malization scale, i.e. the scale at which one performs the measurements). 
The two equations are formally identical at first order in pertubation the- 
ory]^. This is not surprising because one can always define the renormalized 
coupling at a scale a to be equivalent to the bare coupling with a cut-off 
a = a. 

In any practical lattice simulation one computes the regularized integrals, 
eq. (^3l), for a finite a, using the Monte Carlo technique. The coupling 
constant is an input parameter that fixes the scale of the problem. Relating 
the values of the coupling constants at different lattice spacings is, in general, 
a fine tuning problem. 

3.3 Improving the convergence 



It has been shown by Symanzik ||2^ that it is possible to improve the con- 
vergence of the correlation functions of a regularized theory to its continuum 
limit (a — >■ 0) from 0{a) to 0{a"'^^). In order to achieve this, it is necessary to 
add to the action S-^ terms which are proportional to aya"^, ...ja"^ and adjust 
the corresponding coefficients. These tersm are called irrelevant operators 

^Only at second order in perturbation theory the RGE for the renorniaUzed parameters 
shows a dependence from their exact definition, the renormaHzation prescription. 
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since their contribution vanishes in the hmit a ^ 0. This improvement tech- 
nique is heavily used in lattice simulations where the minimum length scale, 
the lattice spacing a, cannot be reduced arbitrarily, therefore one desires the 
dependence on a of the correlation functions to be as small as possible. 

3.4 Lattice regular izat ion and momentum cut-ofF 

We have shown as any given smooth field 0(x) defined in [0, L] can be ap- 
proximated with S{a,x) functions 

K-l 

(f){x) ^ 0^''"(a, x) =^ ^ 0fe(5(a, X - ka) (56) 

fc=0 

with K = L/a. 

The same field expanded in Fourier components 



la;) 

n=0 



r-L 



:) = $^&ne'^"" (57) 

where 

Pn =' ^; b^ =' ^ r 0(x)e---dx (58) 

Also the right hand side of eq. (^) can be expanded in Fourier components 
and this can be written as 

oo 

0i^"(a,x) = ^6;e^P"^ (59) 

n=0 

where 



1 ^-1 r pL 

b' =—Y 

k=0 



(pk I S{a, ka - x)e '^"''dx 





(60) 



It becomes evident that for p„ > l/a the integrand oscillates fast and the 
corresponding integral, b'^, has to be small; while for pn < l/a the integral is 
almost constant and approximately equal to e^*^"'^", therefore b'^'::=^bn- The 
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Figure 10: Different behavior of the integrand of 6^ for high frequency modes 
(left) and low frequency modes (right) respectively. The x axis is in units of 
x/L. 



different behavior of the integrand is shown in figure |T0|. This proves that 
eg .(^61) can be written as 



0(x)~0^^"(a,a;) 



b'°{a,x)=^ 



u\ 



n=0 



Pn bnC'P" 



(61) 



If (f){x) describes some physical quantity and one has a finite resolution 
in space, a, then 0(x) can be replaced by its Fourier expansion with a cut- 
off in momentum space Pmax < f^"^. Therefore a momentum cut-off can be 
regarded as an alternative procedure to regularize the path integral. 

Note that bo is the mean value of 0(x). Moreover pi = 2tt/ L represents 
the minimum energy /momentum mode that can propagate on a finite unidi- 
mensional volume of length L. 

The superscripts "latt" and "co" , used to identify the two different regu- 
larizations, are abbreviations of lattice and cttt-Oj^ respectively. 



3.5 Remarks on the physics of Effective Theories 

We have seen how, in QFT, one is forced to introduce three different length 
scales: 

• a the higher spatial (temporal) resolution of current experiments. 

• a: the scale that corresponds to the precision of the mathematical 
description, the cut-off in the Kadanoff- Wilson approach to Renormal- 
ization Group. The theory is said to be renormalizable if this scale 
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can be arbitrarily small. This does not mean that one can trust the 
predictions of the theory with arbitrary precision. 

• 1/A: the typical lenght scale of the physics being studied. This scale 
is in nature and there is no freedom to fix it. For QCD it is Aqcd — 
250MeV (1/Aqcd ^ Ifm). 

It is usually possible to model the same phenomenon using different QFTs, 
which differ in the regularization-renormalization prescription and/or the 
renormalization scale. The predictions of these different theories must be 
compatible with each other apart from order 0{a) corrections. 

Some QFTs have a finite number of coupling constants and it is possible 
to give a well defined meaning to the limit a — ;> because all the possible 
divergencies can be absorbed in the renormalized constants. These QFTs are 
said to be renormalizable. Other QFTs are not renormalizable because it is 
not possible to absorb all the divergencies in a finite set of constants. The pos- 
sible correlation functions, at different orders in perturbation theory, exhibit 
an infinite variety of divergent behaviour. Originally it was believed that 
only renormalizable QFT made senseQ. The modern picture is different: if 
the world is described by a continuous QFT, it must be a renormalizable one. 
But the world could have a minimum length scale and the renormalization 
requirement is no longer a fundamental one. Perhaps the most important 
modern interpretation of these results is that, if one wants to formulate a 
QFT to describe physics down to a finite resolution, a, and one does not 
pretend it to be the theory of everything, it does not have to be renormaliz- 
able (because one does not pretend to send the scale a to zero) |18[. These 
particular kind of theories are called effective theorie^. 

In the real world, there might be new supersymmetric interactions or 
superstring, or electrons and muons may have internal stucture, none of 
which is incorporated in the Standard Model. Nevertheless the Standard 
Model, as it is today, explains the results of our experiments with a typical 
accuracy of 10%. Renormalization saves us by saying that one does not 
need to know what happens at short distance (high momenta) in order to 
understand low energy experiments. 



^"This requirement led Weinberg, Glashow and Salam to formulate the Standard Model. 
^^E.g. the Fermi theory of electroweak interactions is not renormalizable but it is able 
to describe with good accuracy weak interactions at energy scales below mw 
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Since we will never be able to probe the physical world at every length 
scale, every quantum field theory should be considered an effective theory. 
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4 Lattice QCD 



He who can properly define and divide 
is to be considered a god 

(Plato) 



At the typical hadronic scale, while electroweak effects can be computed 
in the standard perturbative way, QCD effects are not. Therefore it becomes 
necessary to formulate QCD in such a way to make it possible to perform a 
numerical computation of the correlation functions. 

To reach this goal we need to regularize the theory on a lattice by intro- 
ducing a finite lattice spacing a and write down the QCD action in terms of 
the discretized fields. 

In this section we will present one possible discretization of the action 
(due to Wilson) and the problems connected with the definition of chiral 
(massless) fermionic fields on the lattice. 

In principle one can compute the correlation functions of QCD with arbi- 
trary precision by reducing the lattice spacing a. In practice one only has a 
finite computational power therefore one cannot arbitrarily reduce a. Hence 
we will show how to "correct" the action to improve the correction of the 
correlation functions from order 0{a) to order 0(0^). 

We will then discussed the errors associated with the typical lattice ap- 
proximations. 

In the end we will present, as an example, a full program to compute 
Ib^^itT'B and some recent lattice results for this quantity. 

4.1 Basic degrees of freedom and action 

The Aharonov-Bohm experiment revealed that the gauge field A^ = t'^A'l^ is 
not observable, because it is gauge dependent, but the phase of a particle in 
a gauge field background, moving from x to y, is an observable 

-pe^aF^A^d^^ (62) 

{V indicates a path-ordered exponential). 
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For a finite lattice spacing a, it becomes convenient to write the action of 
QCD in terms of the shortest paths on the lattice 

U^{x) =^ e''^^'"'" ''-''''' -l+tgaA^{x + aJl/2) (63) 

t/_^(x) =^ ei9/:-'"'A^dxM _ f;t(a; _ a^) ^ 1 _ igaA^ix - a/2/2) (64) 

(which are associated to the straight path connecting two consecutive lattice 
sites). In QCD U are SU{3) matrices, called links. 

The fermionic degrees of freedom instead, the quarks, can naively be 
defined as q^ix) fields associated to the lattice sites x (where a is the spin 
index and i is the color index). 

The basic discretized operators which appear in the Lagrangian, are: 

• Ordinary derivative^: 

d,j,q{x) = — [q{x + ajl) - q{x - 041)] (65) 

• Covariant derivative: 

def 1 

D^tq{x) = — [U^{x)q{x + a]2) - U-^{x)q{x - op)] (66) 

It is trivial to check that in the continuum limit it is equivalent to the 
usual covariant derivative of QCD. In fact, up to order a corrections, 

D^,q{,x) = —[{l + igaA^{x) + ...){l + ad^ + ...)q{x)- 

(1 - igaA^{x) + ...)(1 - ad^ + ...)q{x)\ 
= {d^ + igA^{x) + ...)q{x) (67) 

For practical purposes, that will be examined below, one is usually in- 
terested in the Euclidean formulation of Lattice QCD. This is achieved by 
performing the Wick rotation]^ 

xq -^ ixo; Xi -> Xi (68) 



^^We remark that a different choice could be made as long as in the limit a ^ one 
recovers the ordinary derivative. 
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The way how different quantities transform under the Wick rotation is reported in 



table ( llOlD 
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Under this rotation the exponential term in the action becomes real 

^iS ^ ^ia^ E. ^ _^ e-5B _ g-a^ E. Ce (gg) 

From now on all the quoted quantities (including 7 matrices) will be Eu- 
clidean. In terms of the links, any correlation function in QCD can be written 

as 

(0|C(...)|0)'""''= /"[df/][dg][dg]C(...)e-'5^[^''''^"l (70) 

where [^ 

5e = 5f "^' + S^''^''' (71) 

and the two contributions are the gauge (pure Yang-Mills) and quark respec- 
tively. We write them in terms of the links as 



,4 



1 - -Re trP^,(x) 



J2h'^%i^^G^'"'(^) + 0ia') (72) 



with 



and 



Pfiu = U^{x)U^{x + afi)U^^{x + aft + ai))U_t,{x + ai>) (73) 



X 

The variable (3 = 6/g'^{a) has been introduced to conform to the standard 
notation of Lattice QCD. In the gauge part of the action there are no order a 
corrections and the first corrections arise at the order a^. On the other side, 
in the quark part of the Lagrangian, order a corrections play, in general, a 
very important role for the following two reasons: 

• If one neglects order a corrections and naively uses the lattice covariant 
derivative to implement iS^"^'^ , one obtains a free quark propagator of 
the form 

S{P) = . , . / ,^ (75) 

«7^ sin[p^a) + am 
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which has 16 zeros in the Brillouin zone in the hmit m ^ 0, to be 
confronted with the single zero of the continuum propagator 

Sip) = -. ^ (76) 

This problem is known as doubling. To get rid of this proliferation of 
zero modes Wilson proposed to add to the action a term of order a of 
the form 

-a^- ^ g(a;)— [f/^(x)g(a; + a/i) - 2q{x) + U-^{x)q{x - aft)] (77) 

X,fJ, 

In practical simulations r is fixed to be 1, nevertheless it is convenient 
to show its dependence. Note that eq. ( |77D is a mass term therefore it 
explicitly breaks chiral symmery. 

• When computing correlation functions one is always interested in the 
limit a — i> and one would like to improve the convergence of correla- 
tion functions from order a to order a^. This can be done by adding 
to the Action terms of order a that compensate for the discretization 
errors up to the same order |]2^. In particular one can choose a term 
of the form 

-«'^ E Qi^h'^G^A^U^) (78) 

X,fl>U 

where G^,y is a discretized version of the chromo-electro-magntic tensor 
and the constant csw must be fixed somehow. This term is usually 
referred to as clover term. It is important to observe that csw is not 
a new free parameter of the theory, it is uniquely determined by the 
value of P (i.e. by the lattice spacing). 

Including these 0{a) corrections, the quark part of the action can be 
re-written as 

sr^'iu, q, ?1 = ^ E <ii^)Q^y[uUy) (79) 

y 
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where 

Qxy[U] = 6^y -K ^ [{r - Y)Uf^{x)6^,y-f, + (r + 7'')t/_^(x)(5^,y+^] 

J2rYG^,{x)5,y (80) 



KCsw 



2 



/lyu 



and, at tree-leveQ 



/3 = ^; ^ = ^ ^^ . ; csw^ = i (82) 

g'^ya) 2ma + 8r 



The action (|79| - |80D is referred to as Sheikoleslami-Wolhert action (SW) p3| . 

The coefficient a? /2k, can be absorbed in the definition of the fermionic 
fields, q and q. 

In practice, in lattice simulations, a is not an input parameter, since it 
is dimensionfuU. Therefore one fixes the lattice spacing, using dimensional 
transmutation, by setting a value for (3 = 6/g'^{a). One procedes in the 
following way: 

One measures on the lattice some dimensionfuU quantity, m (for example 
the K mass or the charmonium IP — IS" splitting), in adimensional units am 
and compares it with the experimental value m'^^^'. From the comparison 
one determines the value of a corresponding to the arbitrary input value of 
p. So that one adjust a by fine tuning the bare /3. 

It is not surprising that the only way to fix the physical parameters of 
the theory is by comparing them with physical experiments. 

On the other size, the improvement coefficients, csw for example, and 
Kcrit (the value of k associated to chiral fermions) are uniquely associated 
with the value of (3 and can be determined theoretically. 



^''Thc fact that Lattice QCD with the action of eq. (JS^) is 0{a) improved for on-shell 
quantities does not simply appear form a Taylor expansion. To show the improvement it is 
necessary to list all dimension 5 operators, use the equations of motion to reduce some of 
them and absorb the contribution of the remaining two operators in the coupling constant 
and mass renormalization 

g ^g {1 + bgma); m —^ m{l + bmma) (81) 
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4.2 Simulation aspects and quenching 

First of all, one can consider correlation functions that do not depend on 
quark fields and integrate out the quarks 

(0|O(...)|0)i"" = /"[dt/][dg][dg]0(...[f/])e-'^^[^''''5"l 

[dU]0{. ..[U]) det Q [U] e^-^^'"'" '^l 

[dU]0{...[U])P[U] (83) 



where 

P[U] = e-'5r"'[C^]+lndetQ[C/] ^g^^ 

In practice one neglects the contribution of Indet Q[?7] in the probability 
P[f/], eq. (0). This is called the quenched approximation. This is a very 
crude approximation and it breaks the unitarity of the theory. Its only moti- 
vation is the limitation in present computer power. It introduces a systematic 
error in the computations that has to be quantified. 

The probability P[U] is real (this is why Lattice computations are per- 
formed in Euclidean space) and resembles a Boltzman weight factor. There- 
fore standard statistical mechanics techniques can be applied. 

We have seen how, for a finite lattice spacing, the Path Integral reduces 
to a well defined i^-dimensional integral which can be computed numerically 
using Monte Carlo integration. 

In a typical Lattice QCD simulation one discretizes the space-time on 
grid of about 48 sites in the temporal direction and 24 sites in each spatial 
direction. Therefore the computation of each correlation function correspond 
to the computation of a i^-dimensional integral with K = A8x 24^ = 663552. 
This is quite a big integral! 

Any standard lattice simulation begins with the creation of an ensemble 
of gauge configurations {[/t*^}. It is created through a Markov process 0, 
i.e. each configuration f/W is generated from the preceding one, t/['~^l, using 
a Monte Carlo algorithm satisfying the condition 

P(t/[^-i] ^ t/W)P[f/l'^-^l] = P(f/[^1 -> [/I*-il)P[t/W] (85) 
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State 


jG 


jPC 


Operator J 


scalar 


1 — 


0++ 


qq' 




1 — 


0++ 


qjy 


pseudoscalar 


1 — 


0-+ 


q-f^q' 




1 — 


0-+ 


q^o^^q' 


vector 


1+ 


1-- 


qYq' 




1+ 


1~ 


qYl\' 


axial 


1 — 


1++ 


qYl^q' 


tensor 


1+ 


1+- 


g-yM^ig' 


octet 




1- 

2 


(g^^7'7°g'^)(7^r)^..fc 




2 


1- 
2 , 


{q^'l'l'lSW)e^,k 


decuplet 


3 

2 


3 + 
2 


{q^'l'lHq''){q"')e^,u 



Table 1: Example of currents used on lattice and their relative quantum 
numbers, g, q' and q" are different flavours. The superscripts i, j and k are 
color labels. 

where P{U -^ U') is the probability of generating the configurations U' from 
the configuration U . An example of such an algorithm is the Metropolis 
Algorithm. A more sofisticate one is the heatbath algorithm. 

Note that P[U] depends on /? = 6/g'^{a) which is the parameter that fixes 
the lattice scale. The initial configuration f/M is usually chosen to be "cold", 
i.e. when all its links are the identity, or "hot" , when each link is a random 



SU{3) matrix fTT], p5| . 



The computation of any correlation function, as defined in eq. (|5BD, can 



be approximated, in analogy with eq. (^), as an average over the ensemble 
of gauge configurations 

(0|O(...)|0r*:^i5^O(...[f/[^]]) (86) 

i 

where N is the number of generated configurations. 

4.3 Correlation functions and fermions 

The typical quantities that are measured on the lattice are the two and three- 
point correlation functions between currents and their (discrete) Fourier 
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transforms at zero momentum, eqs.(^-|T^) PB| , 

C2{t.) = 5^(0|J°(x)J°t(o)|o) (87) 

X 

Csoit.,ty) = 5^(0|J°(x)O(0)J°t(-y)|0) (88) 



x,y 



Since the lattice metric is Euclidean, the asymptotic behaviour of the 
spatial Fourier transform of the two point correlation function, in the limit 
t^ — i> oo, is given by 

C2{t^) ~ Z^e"'"^*^ (89) 

tx—*00 

where mj is the mass of the lightest state |lj) created by the current J^ and 

. mj'{o)\ij)\ 



V^rrij 



(90) 



From the measurement of C2{t^) and its fit to (PUJ), it is possible to extract 
masses of particles, mj. In the same fashion from the asymptotic behaviour of 
the ratio between the three and two-point correlation functions it is possible 
to extract matrix elements ET 



C,0{t.,ty) ^ 1 {lj\0\lj) ^g^^ 



r-^ 



C2{tj)C2{tj,) t^,ty^^ Zj 2mj 

The most general current J^{x) is expressed in terms of fundamental 
fermionic fields ql,{x) (the quark fields). A list of some interesting currents 
J is reported in table ^ In expressions like eq.(^) and ([88| ) these fields are 
Wick contracted 

sz{x,y)=Umii^),i{y)m (92) 

Despite the fact that fermions are neglected when gauge configurations are 
created, they are re-introduced at a later stage as particles propagating in 
the gluonic background field. Therefore the two and three point correlation 
functions can be written as appropriate traces of propagators, S^^^{x,y, [L'"]), 
in the backgroud gluonic field U . For example the propagator of a heavy-light 
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pseudoscalar meson (associated to the current J^ = /i*7'^7^g*, where h, q are 
the heavy and the hght quark respectively), from to x can be computed as 

(o|j°(x)j°t(o)|o) = {om^){i'i'rq:ix)qimi'i'f'Kmo) m 
= {o\su^MiVrsi\M^)iivr\o) 



{U} 



%i^, 0, [u]){i'i'rsLM ^, [u]){ivr 



(i,j are color indices and a, b, c, d are spin indices). Then by making use of 
the H discrete symmetry (see Appendix A) and properties of the 7 matrices 
one obtains 

{Q\j\x)j'\m) = ^Y.^^{sU^Mu])sL,{xMu])] (94) 

{u} 

(the trace is only on the color indices i and j). 

On each gauge configuration f/, the fermion propagator S is computed 
by inverting numericaly the fermionic matrix 

S{x,y,[U]) = {Q[U])-l (95) 

This is the most expensive part of any lattice calculations. In the computa- 
tion of the propagator k and csw are input parameters. The former is in one 
to one correspondence with the fermion mass (the pole-mass) 

m = -hi(l + -(-~^\\ (96) 

a \ 2 \k KcritJ J 

and Kcrit is a parameter depending on /3. The chiral limit corresponds to 
the limit k, -^ Hcru, when the quark becames massless. In practice any 
inversion algorithm for eq. (|95|) converges slower and slower as the chiral limit 
is approached and this can never be reached. 

There are two standard ways of computing the fermion propagator: exact 
and stochastic. The former is very time expensive therefore it is usual normal 
practice to compute propagators ending in one single point of the lattice. The 
latter is less precise but allows one to compute propagators from each point 
to each point of the lattice in a feasible time. 
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4.4 Lattice errors 

Numerical simulations of Lattice QCD are characterized by a number of 
statistical and systematical errors which will have to be taken into account 
when quoting lattice results. What follows is a list of the most common 
errors one has to consider, possibly reduce and, hopefully, quantify: 

• Discretization errors a. Physical results are extracted from lattice 
in the limit of a -^ 0. This limit cannot be reached in real lattice 
simulation and in practice one performs simulations with a finite (as 
small as possible) lattice spacing. The discrepancy between the com- 
puted correlation functions and their continuum limit is usually of the 
order of a (or a^ for improved lattice actions). In typical simulations 
1/a = 1 ^ 3GeV. 

• Statistical errors. All Monte Carlo simulations are based on statisti- 
cal sampling therefore they introduce a statistical error that is expected 
to decrease with 1/yN where N is the number of independent mea- 
surements (in case of one measurement for gauge configuration, A^ is 
the number of uncorrected gauge configurations). Since the gauge con- 
figurations are created using a Markov chain based on small changes 
to link variables, one may be worried about correlations among the 
different configurations. Because of modern day computational power 
it is possible to generate reasonable statistical samples and this is no 
more a major problem. 

• Finite volume. Because of the finite volume, periodic or anti-periodic 
boundary conditions are imposed for the field. Therefore every observ- 
able which is computed on the lattice suffers from an unphysical contri- 
bution of mirror states. In any case, these finite-volume contribution to 
the correlation functions falls off exponentially with the lattice length 
and they are usually negligible for a lattice size L bigger than b/niT^. 
On the other side the effects of mirror states are crucial in preventing 
a direct determination of scattering phases from lattice 



Quenching. This approximation is the hardest to justify. Its only 
reason is the present limitation in computational power. As a conso- 
lation one can argue that present exploratory unquenched simulations 
suggest that the effects of the quark loops in the mass spectrum are 
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small, but the error introduced by quenching in the determination of 
a~^ can be as big as 10%-20%. 

Chiral extrapolation. It has been shown that because of the finite 
volume effects, an infrared cut-off is naturally associated with the lat- 
tice, therefore the u and d quarks are too light to be simulated, even 
by modern day computers. Therefore one usually performs lattice sim- 
ulations for values of rriu = rrid much bigger of the physical values 
(typically bigger than 50 MeV), then performs an extrapolation of the 
results to the limit rriu = rriii = 0. This extrapolation is called the chi- 
ral extrapolation. It corresponds to the limit k, —>■ Hcrit For the masses 
of light particles (the pseudo-Goldstone boson) this extrapolation is 
guided by predictions of the Chiral Lagrangian such as the Gell-Mann- 
Okubo formula. 

Heavy quarks. The c and b quarks are very heavy therefore not all of 
their modes can propagate on a typical lattice. To solve the problem 
there are three common approaches. One possibility is to simulate 
these quarks with a mass smaller than the physical one and then to 
perform an extrapolation to the physical mass (guided by the Heavy 
Quark Effective Theory). The second possibility is to implement the 
HQET on lattice. This implies that one considers the heavy quark 
as static (non-relativistic) and, in principle, systematically computes 
corrections to this approximation in the l/nih expansion. The third 



approach is also based on the HQET and is explained in ref. 28 



Matching between lattice and continuum scheme. Experimen- 
tal data are analyzed using some continuum renormalization scheme, 
usually dimensional regularization with the MS prescription. To con- 
front Lattice QCD results with phenomenology it is therefore necessary 
to match the matrix elements between the two different schemes. In 
general 

(0|O,(...)|0)^ = Zi^{0\Oj{...)\oy^'' (97) 

= (5,, + O(a,))(0|O,(...)|0)'^** (98) 

where Zij are called matching coefficients and have a perturbative ori- 
gin. The matching coefficients can be computed in perturbation theory 
and usually they are known only at 1-loop. Since a^ is big at the typical 
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lattice energy scale, corrections of higher order in a^ can contribute to 
an error in the matching of as much as 10%. Moreover in the match- 
ing procedure it is common that matrix elements of some continuum 
operator mix with the corresponding matrix elements of new opera- 
tors that appear on lattice, because Zij is, in general, non diagonal. 
The contribution of these operators can be big and must be taken into 
account. 

4.5 One full application, /s^/me 

We have seen in section 1 how one could extract fsy/inB by computing 
the Path Integral on the right-hand side of eq. (^, the 2-point correlation 
function, and fitting the result with 

C,{t) = IkHRe-"^^' (99) 

Program C2.C in the Appendix is an example of a real QCD program 
(written in C++) that computes C2{t). It is based on the mathematical library 
described in ref. pOf . 

This program is organized in the following way: It opens the libraries and 
declares the variable containing the parameters of the lattice simulation: 

• beta= /3, that fixes the size of the lattice spacing; 

• mq= m,, the pole mass of the light quark; 

• mh= TUh, the pole mass of the heavy quark; 

• grid_size, the lattice size (for example 16 x 6'^); 

• lattice, the object that contains the grid and its properties (including 
the functions to move on the lattice, a local random number generator 
and functions necessary for the parallelization). 

It then defines the basic fields: 

• U(x,mu) (ij)= UJj'{x), the gauge field configuration; 

• Sq(x,a,b) (ij)= (0|{g^(x), g^(0)}|0), the light quark propagator; 

• Sh(x,a,b) (ij)= (0|{/i^(a;), /i;^(0)}|0), the heavy quark propagator; 

42 



(they are fields of 3 x 3 color matrices and a and h are spin indices) . 

The program starts with a hot configuration (set_hot(U)) then computes 
and discards the first 100 gauge configurations of the Markov chain. The 
Heatbath (heatbath(U)) algorithm is used to generate the Markov chain. It 
is equivalent but more sofisticate than the Metropolis algorithm. 

Then, each 10 gauge configurations of the chain, the program computes 
the light and heavy propagator and measures the 2-point correlation function 
C2(t), eqs. (I)® and (0), 

5^J°(a;)J°t(0) = 5^tr{5,„,(x,0,[[/])5L,(x,0,[[/])} (100) 

X X 

(t = Xo is kept fixed in the sum). 

This is done in the following piece of code: 

f orallsites(x) . . . 
t=x(0); 

F(conf ig,t)+=real(trace(Sq(x,a,b)* 

hermit ian(Sh(x, a, b)))) ; 

Note that the average of the trace in question is always real by virtue of 
the theorem of Appendix A (applied to the case p = 0). This is in agreement 
with naive expectations from eq. (|). 

Finally the program computes the average over the gauge configurations 
of C2(t) (with its the Bootstrap error) and prints out the results. 

The output of the program is plotted in fig. ^ for some arbitrary input 
values of the parameters. In this example, in fact, we did not attempt to tune 
the parameters properly since our main concern was to have a fast running 
program for didactic purposes. 



For a real state-of-the-art computation of fB^/^^^^ we refer now to ref. ||29 
The program used in that simulation is very similar to the one discussed here. 
The only operative differences are in the choice of the parameters and in the 
size of the lattice. Fig. [TT] shows the results from ref. [^ 

The values of /sy^me are computed by fitting C2{t) for different sets of 
input parameters niq and /?, while ruh is tuned to the B quark pole-mass. 
These results are first extrapolated to m^ — ^ and then extrapolated to 
a^O {P ^ oo). 

The results obtained in this way are renormalized in the lattice scheme 
at the lattice energy scale a~^. To obtain the numbers renormalized in the 
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MS scheme at the m^ energy scale they need to be corrected by a matching 
factor. 

In the paper in exam the lattice spacing a is measured (as function of 
(3) by confronting the numerical result for the IP- IS mass splitting in the 
charmonium spectrum with experimental results. 
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Figure 11: Determination of fB^/rn^, fD^/rrlB,fBs^/>^ and /d» y/mn, (Fer- 
milab [^). The two point correlation functions, eq. (D are measured for 
different values of the lattice spacing a and different masses for the light 
quarks m = l/aln(l + 1/(2k) — l/{2nc)) and fitted to extract fsy/^TT-B- The 
results are extrapolated to the chiral limit m —* (top-left). Then the lattice 
spacing is determined by measuring the IP-IS mass splitting in the charmo- 
nium spectrum (top-right). The chirally extrapolated values for /By^me, 
fD^/nfiDjEs \/^nfiBs ^^d /ds ^/^nfiri^ are extrapolated to the continuum limit 
a -^- (bottom- left and bottom-right). 
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A Euclidean Space-Time in o? = 4 dimension 



A.l Wick rotation 

The Euclidean action is obtained from the Minkowskian one by performing a 
Wick rotation. Under this rotation the basic vectors of the theory transform 
according with the following table (E for Euclidean, M for Minkowski) 



E 


M 


E 


M 


x' 


2X0 


x' 


x' 


go 


-ido 


Qi 


d^ 


A' 


-IAq 


A' 


A, 


FH 


-^Foi 


pij 


F^, 


7° 


7° 


i' 


—27* 


7^ 


l' 







and the integration measure transforms as follow 



exp(— 5^;) = exp{iS 



M 



where 



Sf 



d xeCeI--] 



/ d'^XMjCMl--] 



The choice d'^XE = id'^XM can be made, hence ££;[...] 
The Euclidean metric tensor is defined as 



-^m[- 



9e 



-<5^^ = diag(-l, -1,-1,-1) 



(101) 



(102) 



(103) 



(104) 



A. 2 Spin matrices 

• Dirac matrices (Dirac representation) 



7 



1 
-1 



r 



—i(7i 
iai 



7 



1 

1 



(105) 



Dirac matrices (Chiral representation) 

—ia,: 



7 



1 

1 )' ^^ 



iaj 



7 



-1 
1 



(106) 
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All the Euclidean Dirac matrices are hermitian. The following relations 
hold 

(107) 

(108) 
(109) 









9^" = 


2^^ 


,7^} 


= 51^^ 








a"^ = 


2^7 , 


7l 










l' = 


7V7V 




and all the 


(jt^y 


are 


hermitian. 








Projectors 




















L-'- 


1' 


R = 


1 + 7^ 



;iio) 



Traces 



t^i^Y) 


= 4(5^^^ 




(111) 


tr(7'^7'^70 


= 




(112) 


iT{YYYY) 


= A{5^"'5P'' - 


- b^Pb""" + 5^" 5"^) 


(113) 


tr(7^7'^7^7V) 


= A^j^P'' 




(114) 


where e^^^^ = — 1. 









A. 3 Lattice discrete symmetries 

The lattice formulation of QCD is invariant under the following discrete 



symmetries of the quark propagator ||3TI 

• Parity, P: 

S%{^.V. m = l%S%p.{x'',y^, [U''])^%, (115) 

• Charge conjugation, C: 

SU^,y, [U]) = hVU'Si:,,{y,x, [U^]){lVhfs (116) 

• Time reversal, T: 

SU^,y, [U]) = {l'l%a'S%{x^,y\ [U^]){l'l\', (117) 
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• H symmetry: 

S%{x,y. [U]) = lk'Si),,{y,x, [U]hl,^ (118) 

U^, U'^ , 11"^ are the parity reversed, charge conjugate, time reversed gauge 
configurations respectively. 

A. 4 Theorems about correlation functions 

These discrete symmetries play a very important role because they put con- 
traints on the Euclidean correlation functions. In particular one can consider 
correlation function of the form 

G'("''")(pi,...,p;) = /"dVe'Pi"i...dVe'^""" (0|tr{^7^^^7^^..^7'^™} |0) 

(119) 

where the trace is in spin and color, and S are quark propagators connecting 
an arbitrary couple of points in the ensemble {xi, ...,Xn}- Imposing invari- 
ance under P (parity), one obtains that 

G("''")(pi, ...,p-;) = (-l)^G("'"^)(-pi, ..., -pn) (120) 



where A^ is the number of indices /ii,...,/im that differ from 0. Eq. (|120D is 
true also in the Minkowski space. 

Imposing invariance under PCH one obtains that 

G("''")(pi,...,p-;) = (-1)^ [G("''")(pi,...,p-;)]* (121) 

This tells whether any correlation function is real or imaginary. 

Eq. ( |121| ) is not true in Minkowski space. It is replaced by an equivalent 



expression where N counts the total number of indices /ii,...,/im that are 
equal to 5. 

As an example we consider 

G(3.2)(p^ = / ^i^^ip^ (o| tr {S{x, 0)7^7^5(0, x)-f^} |0) (122) 



Since A^ = 3, eq. ( p.20|) tells that it is odd under parity and eq. ( |121| ) tells 



that it is imaginary. Hence for p = it must be zero (because it is odd under 
parity). 
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B Example programs 

B.l Basic integration: programl.c 

// programl.c 
#include "random. c" 



double f (double x) { 
return sin(Pi*x) ; 

>; 

int mainO { 

double Sum, 11, 12, 13, xO, xl, x2; 
double alpha=0 ; 
double beta=l; 
long N=20,i; 

// ////////////////////// 

// Basic method 

// ////////////////////// 

for(N=2; N<10000; N*=2) { 

Sum=0 ; 

for(i=0; KN-1; i++) { 

xO=alpha+ (beta-alpha) *i/N ; 

Sum=Sum+f (xO) ; 

>; 

Il=Sum* (beta-alpha) /N; 

// ////////////////////// 

// Newton-Cotes method 

// ////////////////////// 

Sum=0; 

for(i=0; KN-1; i++) { 

xO=alpha+ (beta-alpha) *i/N ; 

xl=xO+ (beta-alpha) /N ; 

Sum=Sum+ (f (xO) +f (xl) ) /2 . ; 

>; 

I2=Sum* (beta-alpha) /N; 

// ////////////////////// 

// Simpson's method 

// ////////////////////// 

Sum=0 ; 

for(i=0; i<N-2; i++) { 
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xO=alpha+ (beta-alpha) *i/N ; 

xl=xO+ (beta-alpha) /N ; 

x2=xl+ (beta-alpha) /N ; 

Sum=Sum+ (f (xO) +4 . 0*f (xl)+f (x2) ) /6 . ; 

>; 

I3=Suiii* (beta-alpha) /N ; 

printf ("7.i\t7.f\ty.f\t7.f\n",N, 11, 12, 13); 

}; 

>; 



B.2 Monte Carlo integration: prograin2.c 

// program2.c 
#include "random. c" 

double f (double x) { 
return sin(Pi*x) ; 

>; 

int mainO { 

double 1, Sum,x; 
long N ; 
Sum=0 ; 

for(N=l; ; N++) { 
x=Random() ; 
Sum=Sum+f (x) ; 
l=Sum/N; 
printf ("N = y.i, 1(N) = 7.f\n",N,l); 

>; 
}; 



B.3 Monte Carlo integration in 3D: programs. c 

// programs. c 
#include "random. c" 

double f (double *x) { 

return 3 . 0*x [0] *x [0] *x [1] +2 . 0*x [2] *x [2] *x [2] ; 

>; 

int mainO { 

double 1, Sum, x[3]; 
long N; 



50 



Sum=0; 

for(N=l; ; N++) { 

x[0]=Random() 

x[l]=Rajidom() 

x[2]=Random() 

Sum=Sum+f (x) ; 

I=Sum/N; 

printfC'N = 7.i, I(N) = 7.f\n",N,I); 

>; 
}; 

B.4 Metropolis Monte Carlo integration: prograin4.c 

// progrcLm4.c 
#include "random. c" 

double P (double *x) { 

return exp(-(x[0] *x [0] +x [1] *x[l] +x[2] *x[2] ) ) ; 

}; 

double f (double *x) { 

return 128 . 0*x [0] *x [0] *x [0] *x [1] *x [1] *x [2] ; 

>; 

int mainO { 

double 1, Sum, x [3] , y[3],Py,Px; 
long N.j.step, Nstep=1000; 

Sum=0; 

X [0] =Random() ; 
x[l]=Random() ; 
x[2]=Random() ; 

for(N=l; ; N++) { 

for(step=0; step<Nstep; step++) { 
y[0]=Random() ; 
y [l]=Random() ; 
y[2]=Random() ; 
if(P(y)/P(x) > RandomO) { 

X [0] =y [0] ; X [1] =y [1] ; x [2] =y [2] ; 

>; 
>; 

Sum=Sum+f (x) ; 
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I=Sum/N; 

printfC'N = 7.i, I(N) = 7.f\n",N,I); 

}; 
>; 



B.5 Metropolis and Bootstrap: programs. c 

// programs. c 
#include "random. c" 

long D=3; // number of dimensions 

long N=8192; // number of Montecarlo samples 

long M=100; // number of Bootstrap samples 

long Nstep=1000; // number of interations per config. 

double P (double *x) { 

return exp(-(x [0] *x [0] +x[l] *x[l] +x[2] *x[2] ) ) ; 

>; 

double f (double *x) { 

return 128 . 0*x [0] *x [0] *x [0] *x [1] *x [1] *x [2] ; 

>; 

void swap (double &a, double &b) { 
double c; 
c=a; a=b ; b=c; 

>; 

int mainO { 
double 1, Sum; 
double X [N] [D] ; 
double y [D] ; 
double Ibar [M] ; 
long i,j,d,k[N] [M] , step; 

// ///////////////////////////////////////// 
// generate conf if gurations using Metropolis 
// ///////////////////////////////////////// 

for(d=0; d<D; d++) x[0] [d] =Random() ; 
for(i=l; i<N; i++) { 

for(d=0; d<D; d++) x[i] [d] =x[i-l] [d] ; 
for(step=0; step<Nstep; step++) { 
for(d=0; d<D; d++) y [d] =Random() ; 
if (P(y)/P(x[i]) > RandomO) 
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for(d=0; d<D; d++) x [i] [d] =y [d] ; 

}; 
>; 

// iiiiiiiiiiiiiiimmmiiiiimiiiiiiim 

II compute the numerical integral I(N) 

// llllllllllllllllllllllllllllllllllllllllll 

Sum=0; 

for(i=0; i<N; i++)Sum=Sum+f (x[i] ) ; 

I=Sum/N; 

printfC'N = y.i, I(N) = 7,f\t", N, I); 

// llllllllllllllllllllllllllllllllllllllllll 

II create Bootstrap configurations 

// llllllllllllllllllllllllllllllllllllllllll 

for(i=0; i<N; i++) 
for(j=0; j<M; j++) 

k[i] [j] = (long) ((double) N*Random()); 

// llllllllllllllllllllllllllllllllllllllllll 

II create Bootstrap integrals Ibar[j] 

// llllllllllllllllllllllllllllllllllllllllll 

for(j=0; j<M; j++) { 
Sum=0 ; 

for(i=0; i<N; i++) Sum=Sum+f (x[k[i] [j]] ) ; 
Ibar[j]=Sum/N; 

>; 

// llllllllllllllllllllllllllllllllllllllllll 

II sort the Ibar[j] 

// llllllllllllllllllllllllllllllllllllllllll 

for(j=0; j<M-l; j++) 
for(i=j+l; i<M; i++) 
if (Ibar[j]>Ibar[i]) 

swapdbar [i] , Ibar [j] ) ; 

// llllllllllllllllllllllllllllllllllllllllll 

II print the confidence interval 

// llllllllllllllllllllllllllllllllllllllllll 

printfC and It < I < y.f \n" , Ibar [16], Ibar [84]); 
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>; 



B.6 Marsaglia random number generator: random. c 

/* random. c */ 

// immmmmiimimmmimmmm 

II Remastered version of the generator used by 

// the UKQCD Collaboration. 

// It is based on the Marsaglia generator 

// //////////////////////////////////////////// 

#include <stdio.h> 

#include <math.h> 

#include <complex.h> 

#include <sys/time.h> 

#include <unistd.h> 

#define and && 

#def ine or | | 

#define Pi 3.14159265359 

#define PRECISION le-16 

#define TRUE 1 

#define FALSE 

float Randomdong ijkl=0) { 
static float u[98]; 
static float c; 
static float cd; 
static float cm; 
static int ui; 
static int uj ; 
static int first=0; 
int i, j, k, 1, ij, kl; 

if ((first==0) II (ijkl!=0)) { 

printf ("initializing the Random generator\n") ; 
if( (ijkl < 0) I I (ijkl > 900000000) ) 

exit(l) ; 
ij = ijkl/30082; 
kl = ijkl - (30082 * ij); 
i = ((ij/177) •/. 177) + 2; 
j = (ij 7, 177) + 2; 
k = ((kl/169) % 178) + 1; 
1 = kl 7. 169; 
if( (i <= 0) II (i > 178) ) 

exit(l) ; 
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if( (j <= 0) II (j > 178) ) 

exit(l) ; 
if( (k <= 0) II (k > 178) ) 

exit(l) ; 
if( (1 < 0) II (1 > 168) ) 

exit(l) ; 
if (i == 1 && j == 1 && k == 1) 

exit(l) ; 
int ii, jj, m; 
float s, t; 
for (ii = 1; ii <= 97; ii++) { 

s = 0.0; 

t = 0.5; 

for (jj = 1; jj <= 24; jj++) { 
m = ((i*j y. 179) * k) 7. 179; 

i = j; 

j = k; 

k = m; 

1 = (53*1+1) 7. 169; 

if (l*m 7. 64 >= 32) s += t ; 

t *= 0.5; 

>; 

u[ii] = s; 

>; 

c = 362436.0 / 16777216.0; 
cd = 7654321.0 / 16777216.0; 
cm = 16777213.0 / 16777216.0; 
ui = 97; 
uj =33; 



f irst=l; 

}; 

float luni; /* local variable for Float */ 
luni = u[ui] - u[uj]; 
if (luni < 0.0) 

luni += 1.0; 
u[ui] = luni; 
if (— ui == 0) 

ui = 97; 
if (— uj == 0) 

uj = 97; 
if ((c -= cd) < 0.0) 

c += cm; 
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if ((luni -= c) < 0.0) 

luni += 1.0; 
return ((float) luni); 

>; 

B.7 /b and m^. Program: C2.C 

// //////////////////////////////////////////////////////// 

// Program C2.C written by Massimo Di Pierro (§ July 2000 

// //////////////////////////////////////////////////////// 

// WORKING EXAMPLE of a Lattice QCD program to compute the 

// Euclidean propagator of an Heavy-Light Meson, C2(t) 

// To extract m_B and f_B fit output with 

// 

// C2(t) = 1/2 f_B'2 m_B exp(-m_B t) + ... 

// 

// and extrapolate to 

// mq -> (GeV) 

// mh -> mb (GeV) (the b quark pole mass) 

//a -> (GeV-(-l)) 

// //////////////////////////////////////////////////////// 

// #define PARALLEL 

// //////////////////////////////////////////////////////// 

// open the libraries: Matrix Distributed Processing 1.0 

// (for a description read: hep-lat/0004007) 

// //////////////////////////////////////////////////////// 

#include "MDP_Lib2.h" 
#include "MDP_MPl.h" 

// //////////////////////////////////////////////////////// 

// open the FermiQCD libraries 

// //////////////////////////////////////////////////////// 

#include "MDP_Gauge . h" 

#include "MDP_Fermi .h" 

#define GeV 1 

// //////////////////////////////////////////////////////// 

// main program 

// //////////////////////////////////////////////////////// 

int main(int argc, char **argv) { 

mpi . open_wormholes(argc, argv) ; // open communications 
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// ////////////////////////////////////////////////////// 

II declare parsuneters of the simulation 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 

int Nt=16, Nx=6, Ny=6, Nz=6; // lattice size 

int Nc=3; // set colors, SU(Nc) 

float beta=5.7; // set lattice spacing 

float mq=0.2*GeV; // set light quark pole-mass 

float mh=0.7*GeV; // set heavy quark pole-mass 

// ////////////////////////////////////////////////////// 

// additional parameters (they only depend on beta!) 

// ////////////////////////////////////////////////////// 

float a =0.91/GeV; // lattice spacing 

float cSW =1.57; 

float kappa_c=0. 14315; 



// ////////////////////////////////////////////////////// 
// define gamma matrices in auclidean space 
// ////////////////////////////////////////////////////// 
define_base_matrices("FERMILAB") ; 

// ////////////////////////////////////////////////////// 
// define the grid size on which the lattice is defined 
// ////////////////////////////////////////////////////// 
int grid_size[]={Nt,Nx,Ny,Nz}; 

// ////////////////////////////////////////////////////// 
// associate the lattice to the grid 

// ////////////////////////////////////////////////////// 
generic_lattice lattice(4, grid_size) ; 

// ////////////////////////////////////////////////////// 

// define a gauge field 

// U(x,mu) = exp(i g A(x,mu) ) 

// to the sites of the lattice 

// ////////////////////////////////////////////////////// 

gauge_field U(lattice,Nc) ; 

// ////////////////////////////////////////////////////// 

// define a light and a heavy propagator 

// S = <0l q(x) \bar q(0) I 0> 

// to the sites of the lattice 

// ////////////////////////////////////////////////////// 

f ermi_propagator Sq(lattice,Nc) ; 
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f ermi_propagator Sh(lattice,Nc) ; 

// ////////////////////////////////////////////////////// 

II define a variable site 

// to move on the lattice 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 

site x(lattice) ; 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
II define the number of gauge configurations to be used 
// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
int Nconfig=100; 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
II define some more auxiliary variables 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
int i,j,t, config; 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
II define an object for bootstrap error 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
JackBoot C2(Nconfig,Nt) ; 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
II creating initial gauge conf iguratino 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
set_hot(U) ; 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 

II create and skip 100 gauge configuration 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 

U . param . bet a=beta ; 

heatbath(U,100) ; 

// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
II setting the parameters for the light and heavy quarks 
// llllllllllllllllllllllllllllllllllllllllllllllllllllll 
Sq . param . kappa=l . 0/ ( (exp (mq*a) -1 . 0) *2 . 0+1 . 0/kappa_c) ; 
Sh . param . kappa=l . 0/ ( (exp (mh*a) -1 . 0) *2 . 0+1 . 0/kappa_c) ; 
Sq . param . cSW=Sh . param . cSW=cSW ; 
Sq.precision=Sh.precision=le-7; 

// ////////////////////////////////////////////////////// 

// loop over the gauge configurations 

// ////////////////////////////////////////////////////// 
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f or (conf ig=0; conf ig<Nconf ig; coiifig++) { 

heatbath(U, 10) ; // each 10 gauge configurations 

compute_em_f ield(U) ; // compute electromagnetic-field 

generate (Sq,U) ; // compute light propagator 

generate (Sh,U) ; // compute heavy propagator 

// //////////////////////////////////////////////////// 

// compute the pion propagator by 

// wick contracting 

// 

// C_2(t_x) = 

// \sum_{x} \bar h(x) GcunmaS q(x) \bar q(0) GammaS h(0) 

// 

// as function of t = t_x 

// //////////////////////////////////////////////////// 

for(t=0; t<Nt; t++) C2(0,t)=0; 

f orallsites(x) { 

t=x(0); 

for(i=0; i<4; i++) 
for(j=0; j<4; j++) 
C2(config,t)+= 

real(trace(Sq(x,i, j)*hermitian(Sh(x, i, j)) ) ) ; 

>; 

// //////////////////////////////////////////////////// 

// for each t print out C_2(t) with the Bootstrap error 
// //////////////////////////////////////////////////// 
printf ("\nRESULT FOR C2(t) ((S gauge = y.i)\n", conf ig) ; 
printf ( "==================================\n" ) ; 

printf ("t\tC2\t\t (error) \n") ; 

printf ( " ==================================\n" ) ; 

for(t=0; t<Nt; t++) { 

C2.plain(t) ; 

printf ("7,i\t7.e\t7.e\n", t, C2.mean(), C2.b_err()); 

}; 

printf ("==================================\n\n") ; 

f f lush(stdout) ; 

>; 

mpi . close_wormholes() ; // close communications 
return 0; 

>; 
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C status of Lattice QCD 

In this appendix we will briefly report a very small subset of recent Lattice 
QCD results which have direct phenomenological interest and which provide 
an example of the state-of-the-art in Lattice QCD simulations^. 

A huge amount of work has also been dedicated by the lattice community 
to the study of some theoretical and numerical aspects of field theories and 
the properties of their lattice regulated versions. They include the study 
of different possible discretizations for the Dirac operator (^), study of low 
energy eigenvalues of these operators, the confining properties of different 
Yang-Mills theories and chiral symmetry breaking. These studies have given 
some important insights in the understanding of QCD and constitute the 
foundations on which any simulation of phenomenological interest relies on. 

• Chiral symmetry on the Lattice. The light spectrum of QCD is 
dominated by spontaneous chiral symmetry breaking (in fact the pion 
is both a bound state and a Goldstone boson). At a classical level the 
continuum Dirac operator,^, preserves the chiral symmetry, which we 
rewite as the Ginsparg- Wilson relation 

7^|) +^7^ = (123) 

The chiral symmetry is broken, at the quantum level, by the chiral 
anomaly. 

The lattice regularized Dirac operator, eq. (^), explicitely breaks this 
symmetry and therefore does not provide a satisfatory description of 
chiral physics. 

Since today no-one succeded in writing down a discretized versions of 
chiral fermions because of the famous Nielsen-Ninomiya no-go theo- 
rem [^. Recently two solutions have been found to this problem and 
they are both equivalent to modify the chirality condition, eq (|123|) , 
into 

7^_^ +^7^ = ap-i^lp (124) 

where the right-hand side vanishes in the limit a — i> 0. 



^^The papers quoted below are chosen as examples and we do not aim to provide a 
complete list of references on any of the subjects. 
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The two approaches are known as domain wall fermions P^ and Neu- 



berger fermions |^ . They, from a theoretical point of view, are equiv- 



alent ||35|, |36 |. 



The work in this subject has still to be considered exploratory but 
seems very promising. Moreover it provides a theoretical testing ground 
for other well-estabilished lattice discretizations of fermions, such as 
Sheikoleslami-Wolhert (also known as Clover) and Kogut-Susskind (also 
known as Staggered) fermions. 

At present the most significant phenomenological lattice results are 
computed using Clover (and Staggered) fermions. Here "significant" 
means that simulations with Clover (and Staggered) fermions have been 
performed in a wide range of lattice spacings (IGeV < a~^ < 3GeV), 
on relatively large lattices (up to 48 x 24'^ and even bigger) and for many 
different values of the light quark masses. Large statistical Monte Carlo 
samples (of the order of hundreds of configurations) have been created 
and scrutinized. 

Moreover the most important results of phenomenological interest have 
been reproduced indipendently by different international collaborations 
and agree with each other within the statistical errors. 



Light hadronic spectrum. Fig. |12| shows some of the most recent 



lattice results for the mass specrtum of light mesons and baryons, as 
computed by the CP-PACS collaborations. These results are obtained 
using Clover fermions, extrapolated to the chiral limit, in the quenched 
approximation. The error includes the effect of the chiral extrapolation 
but does not include the unknown effect of quenching. 



Fig. |I3| shows a comparison between quenched and unquenched re- 
sults for the mass of light hadrons as function of the lattice spacing. 
Quenched results are easier to compute therefore have a much smaller 
statistical error. At present their error is dominated by the systematic 
one. The comparison with unquenched results can be used to quantify 
this systematic error. 

Glueball spectrum. One sector in which lattice calculations have ab- 
solutely no competitors (because no models are available) is the com- 



putation of the glueball spectrum. Fig. 14 shows some lattice results 
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obtained in ref . |^ . Some candidate resonances have been observed in 
experiments. 

Confinement. Lattice has been the first theoretical tool to give direct 
evidence of the phenomenon of confinement. This is done through the 
computation of the Wilson loop, from which one obtains the chromo- 
electro-magnetic potential between a couple of static quarks-antiquark. 
This potential (both for quenched and unqueched simulations) is shown 
in fig. |15|. For large distances the behavior of this potential, in the 
quenched approximation, grows linearly because no new particles can 
appear in the vacuum to screen the two static particles. The linear- 
ity is explained with the dual-superducting property of the QCD vac- 
uum [^ . It forces the chromoelectric field of the system of two par- 
ticles to be squeezed into a flux tube connecting them: the string. In 
the unquenched (full) theory, instead, this potential only grows up to 
the point when the string contains enough potential energy to create 
in the QCD vacuum a new quark-antiquark pair that breaks the string 
itself. Hence the potential reaches a plateau in correspondence to the 
treshold energy. 

The deviation between the quenched and the unquenched potential 
has not been observed so far because, for numerical reasons, dynamical 
quark masses are still too heavy and lattice volumes too small. 

The computation of the chromo-electromagnetic potential, combined 
with experimental measurements, provides the best present determina- 
tion of Us (fig. |16D, which is the expansion parameter of any perturba- 
tive QCD calculation. 

Matrix elements and decays. Almost all decays of hadrons can 
be parametrized in terms of matrix elements that encode the non- 
perturbative contribution of QCD. For many of these matrix elements 
lattice computations have been able to produce satisfactory results. In 
table we list, as an example, some of those matrix elements. Some 
other results are not conclusive and occasionally very controversial (for 
example attempts to compute e'/e). 

Matrix elements that include a contribution of final state interactions 
(for example B — > tttt) have so far been outside the reach of lattice 
computations because of the Maiani- Testa no-go theorem p 
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Matrix element 


Process 


Parameters 


(0 Irff^qlB) = JBP^ 


B -^ leptons 


\Vtdl Vts 


(0 cYq\D) = foP^ 


D -^ leptons 




{D cTs K) 


D —^ K + leptons 


\Vcs\ 


{K\-sVu\M^) 


K -^ Mu + leptons 


Ks 


{D cTd Md) 


D -^ Mfi + leptons 


\vj 


{B\bTu Mu) 


B ^ Mu + leptons 


\Vub\ 


{B\hTc\D*) 


B ^ D* + leptons 


\V,b 


{B\qTq B*) 


i? ^ TT + leptons 


\Vub\ 


{K\sTqsTqK) oc Bk 


{K-K mixing) 


\Vul Vts 


(DcTqcTqlD) (x Bd 


{D-D mixing) 




{B\hTqhTqB) oc Bb 


{B-B mixing) 


\Vtd , Vts 


{B\bTbB) 


B kinetic and magnetic energy 




{B\hVqqVhB) 


inclusive B decay 




{K^hVqqVh\k^) 


inclusive A^ decay 





Table 2: Examples of matrix elements usually computed on the lattice and 
related processes {Mg stands for the most general {qq') meson, for example 
TT, p or K; r is the most general spin(8)color matrix). The table also shows 
the VcKM matrix elements that are associated to the processes. 
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Lattice, effective field theories and models. Since lattice can be 
used to compute matrix elements and these can be confronted with pre- 
dictions from effective theories or models, it becomes possible to use the 
lattice results to extract their effective parameters. As an example in 
ref. |^2| the effective coupling of the Heavy Meson Chiral Lagrangian is 
measured on the lattice using a numerical computation for the following 
matrix element 

^W = ^ E f {B\A,{^)\B*)dn^ (125) 

where Af^{x) is the axial current and fix is the solid angle associated to 



the 3D vector x. Fig. |I7| shows a comparison between the lattice result 
for E{r) and a prediction of the Chiral Quark Model. The parameters 
of the model have been adjusted to fit the experimental mass spec- 



trum of heavy mesons. See ref. |^ for further details. The agreement 
between the lattice (a first principle simulation) and the model (conse- 
quence of experimental observations and some theoretical assumptions) 
is remarkable. 

Heavy quarks and CKM Matrix. As we have shown in Section 1 
the Cabibbo-Kobayashi-Maskawa matrix elements have to be extracted 
from a comparison between experimental and theoretical predictions. 
The latter, at present, are computed on the lattice with a non-negligible 
uncertainty. Table § includes a list of those processes that mainly 
contribute to the determination of the CKM matrix elements and have 
to be computed using lattice simulations. 

Fig. 0(top) shows the present constraints on the CKM mixing angle 
in the (p, fj) plane (where p and r] are parameters of the Wolfenstein 
parametrization of the CKM matrix) and should be compared with 
fig. |T8|(bottom) which is obtained using the same experimental and 
theoretical input but assuming a possible future uncertainty instead 
of the "real" present one for the lattice parameters. The future un- 
certainty is based on being able to generate 1000 gauge configurations 
with a lattice spacing of a = O.OSfm and with niT^/nip = 0.4 ^1|] (the 
latter constraint measures how well one is approaching the chiral limit). 
This estimate also assumes that systematic error due to quenching are 
under control. 
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The comparison shows how important it is to invest in lattice sim- 
ulations while, contemporary, investing in experimental facilities. In 
fact for many fundamental quantities the theoretical uncertainty is as 
significant (if not more) than the experimental one. 
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Figure 12: Light hadrons spectrum as computed by the CP-PACS collabo- 
ration \M 
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Figure 13: Hadron masses in full QCD as function of the lattice spacing 
compared with quenched results (with Wilson fermions). rriK is used as 
input. Results form the CP-PACS collaboration 
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Figure 14: Mass spectrum in the SU(3) quenched theory as computed in 



ref. |3^. The scale in the right vertical axis is based on the assumption that 
ro = 410MeV. 
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Figure 15: String tension in quenched and unquenched (two dynamical hght 
flavours at k = 0.1575), extracted from ref. |3^ 
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F2 [HERA] 

ev. shapes [HERA] 

QQ -I- LATTICE QCD 
J/^F -I- Y decays 

e+e~[ev. shapes 22 GeV] 

e+e"[ahad] 
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pp ^bbX 
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Figure 16: Comparison among different determinations for as at tlie Z pole 
mass. 
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Figure 17: Comparison of a Lattice QCD matrix element, E{r), and the 
same matrix element evaluated in the Chiral Quark Model. 



71 



-1 -0,S -0-6 -0,4 -0,2 




-1 -0,S -0,6 -0,4 -0,2 




Figure 18: Current (top) and future (bottom) allowed regions for the point 
(p, f/) , the vertex of the unitarity triangle. The three regions correspond to 
5%, 68% and 95% confidence levels. The future estimate is based on the 



assumptions discussed in ref. [^ 
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